%% q1.1 用matlab软件对所研究问题相关的模型其作共线性诊断；
clc
clear
warning off   %为了关掉警告
x = xlsread('data.xlsx','q1');
x1 = x(:,2);   %读取到单产这个因变量这一列
y = log(x1);    %单产因变量全部取对数   对数的底默认为e
x2 = x(:,3:end); %读取后面六个自变量
x = log(x2);    %自变量全部取对数 

regression(y,x) %多元线性回归
%多元线性回归分析结果：
%F值右边的p值是F检验的sig值，通常小于0.05我们认为模型整体显著
%参数估计里面的每个p值，是T检验的p值，通常小于0.05，认为变量显著
gongxianxing(y,x) %共线性检验

%主成分回归步骤：
% 1.进行多元线性回归分析及共线性诊断；
% 2.若存在共线性，进行主成分分析；
% 3.求得主成分得分，确定主成分个数；
% 4.将因变量对保留主成分进行回归分析；
% 5.回代主成分，得到新的回归模型；
% 6.对回归方程给于专业解释；

% 1. 进行多元线性回归的结果 
%    F检验：是否具有显著的线性关系  F值 p值 具有很强的线性关系
%    自变量对因变量的影响均不显著，但观察到了有方差膨胀因子VIF>10，说明自变量间存在多重共线性。

% 2. 进行共线性诊断
%    又观察到条件指数k，可以认为自变量存在中等程度的多重共线性。0-10较小，10-30中等，30-较大

%% q1.2 请利用主成分回归的方法，从而从经济效应的角度得出最佳的棉花补贴方式。
xishu=jiandan(y,x); %简单统计量和相关系数矩阵
xishu1=pca(y,x); %主成分分析
xishu2=zhuchengfen(y,x); %主成分回归

% 可选择前两个主成分替代原来的6个自变量，根据主成分分析的结果得到方程
% 主成分分析的结论
%   第一主成分大约等于这六个投入要素和的一个常数倍，可以简称为总投入主成分；
%   第二主成分的表达式中系数为正的是劳动利润和农药费（看的是特征向量），系数为负的是种子费、化肥费、机械费和灌溉费，
%   其中劳动利润的符号是正的，且系数为0.6722，相比其他投入占有重要优势，可以简称为药物及劳动主成分。
%   标准化因变量对主成分进行多元回归分析。可得到标准化因变量对主成分的回归方程
%   考察表中投入要素补贴具有可行性，最佳的补贴方式是农药补贴。

%% q1小结：
% 主成分回归的方法是将原来的回归自变量变换到主成分，选择其中重要的主成分作为新的自变量，
% 丢弃了影响不大的自变量，实际上达到了降维的目的，
% 然后用原来最小二乘法对选取主成分后的模型参数进行估计，然后再变换回原来的模型求出参数的估计。
% 优点：简化结构、消除变量间的相关性（可以写在文章末尾，模型的优缺点）
% 缺点：回归方程的解释比较复杂
% 我们通常仅将主成分回归作为分析多重共线性问题的一种方法。根据实际分析效果来评价方法的适用性。

%% q2.1 运用第1题的数据，做单产对于6个要素投入的岭回归模型，画出岭迹图。
% clc
% clear
% warning off
% x = xlsread('data.xlsx','q1');
% x1 = x(:,2);
% y = log(x1);
% x2 = x(:,3:end);
% x = log(x2);
linghuigui(y,x) %岭回归

%RIDGE参数要选取 RIDGEVIF都 < 10的数

%% q3.1 试将GDP和居民消费水平这2个因变量对其余6个经济发展指标做偏最小二乘回归。

%偏最小二乘法回归针对的是对多自变量的回归建模中，特别是在观察值数量少且存在多重相关性等问题，
%用偏最小二乘回归分析建模，比对逐个因变量进行多元回归更加有效。
%与主成分回归不同的是：不仅吸取了主成分回归中从自变量提取信息的思路，
%同时还注意了主成分回归中所忽略的自变量对因变量的解释问题。

w = xlsread('data.xlsx','q3');
w1 = w(:,1:2); %因变量
w2 = w(:,3:end); %自变量
xiangguanxishu(w1,w2); %相关系数
%分析相关系数表中数据，可知自变量之间存在较大的多重共线性，
%考虑到因变量个数为2超过1，因此采用偏最小二乘回归分析。
%在写作的时候可以写：与***成正相关，与**成负相关。
extract(w1,w2); %标准化自变量组、因变量组与成分的回归结果分析
jiaochajianyan(w1,w2); %交叉检验，得出偏最小二乘提取因子解释的百分比
%交叉有效性检验
%我们根据MATLAB程序中交叉有效性检验的结果可知，对自变量组和因变量组提取2个成分对达到要求。
%表给出偏最小二乘提取2个成分解释的百分比，可以看出提取2对成分对自变量组的解释比率为72.04%，
%对因变量组的解释比率为74.7%，这说明建模效果还是比较好的。
pianzuixiaoercheng(w1,w2); %偏最小二乘回归     注意：结果中为 y1 接 常数项  下一行类似
yuce(w1,w2);
%标准化变量与成分变量之间的回归方程
%因变量组与自变量组之间的回归方程

%% q3小结
% 1.偏最小二乘回归提供了一种多因变量对多自变量的回归建模方法，特别当变量之间存在高度相关性时，
%   用偏最小二乘回归进行建模，其分析结论更加可靠，整体性更强。
% 2.偏最小二乘回归可以有效地解决变量之间的多重相关性问题，适合在样本容量小于变量个数的情况下进行建模;
%   变量之间多重相关性经常会严重危害参数估计，扩大模型误差，并破坏模型的稳健性。
%   偏最小二乘回归采用对数据信息进行分解和筛选的方式，有效地提取对系统解释性最强的综合变量，
%   剔除多重相关信息和无解释意义信息的干扰，从而克服了变量多重相关性在系统建模中的不良作用.

%% q4.1 序列y1，y2，y3分别表示我国1952年至1988年工业部门、交通运输部门和商业部门的产业指数序列，由于三个部门的产出是相互影响、相互制约，可以考虑建立VAR模型，建议先对数据取对数。
%clc,clear   向量自回归
data= xlsread('data.xlsx','q4');  %导入数据
date=data(:,1);  %第一列为年份
data1=data(:,2:end);
ld=log(data1);%对矩阵data取对数，以消除线性趋势
[n,p]=size(ld); %求矩阵ld的行和列

%画对数序列图
plot(date,ld,'-');
legend('工业部门','交通运输部门','商业部门');

%单位根检验
x=zeros(n,p);  
h=zeros(p,1);
pValue=zeros(p,1);
stat=zeros(p,1);
cValue=zeros(p,1);
for i=1:p
    x(:,i)=ld(:,i);  %ld矩阵的每一列赋予x
end

%%每一个自变量进行ADF平稳性检验
% ADF单位根检验结果，结果见表。从表中可以看出，Tau的p值都较大，
% 表明该统计量不显著不同于0,所以不能拒绝序列有单位根的假设，即该序列非平稳。
for i=1:p
    [h(i,:),pValue(i,:),stat(i,:),cValue(i,:)] = adftest(x(:,i),'model','ARD','lags',2);
end
fprintf('\n');
fprintf('-------------------------ADF检验-------------------------\n');
fprintf('%4s%12s%15s%15s','h','统计量','临界值','p值');
fprintf('\n');
for i = 1:p
    fprintf('%.4f%13.4f%17.4f%17.4f\n',h(i,:),stat(i,:),cValue(i,:),pValue(i,:)); 
end

%%Johansen协整检验
% 再对协整性进行检验，r=2的假设被接收，故3个变量间至少存在2个协整关系。
[h,pValuet,statt,cValue] = jcitest(ld,'model','H1','test','trace');%迹检验
[h,pValuem,statm,cValue] = jcitest(ld,'model','H1','test','maxeig');%最大特征根检验

%%建立模型
% 根据VAR(p)模型最优滞后阶数准则确定滞后阶数为1
ld1=ld(1:n-1,:); %对1-18组数据建立模型，预留第19组数据检验模型的预测效果
Md0 = varm(3,0); %滞后0阶向量自回归模型
EstMd0 = estimate(Md0,[ld1(:,1),ld1(:,2),ld1(:,3)]);%滞后0阶向量自回归模型估计
summarize(EstMd0); %展示滞后0阶向量自回归模型模型结果
Mdl = varm(3,1);    %滞后1阶向量自回归模型
EstMdl = estimate(Mdl,[ld1(:,1),ld1(:,2),ld1(:,3)]);%滞后1阶向量自回归模型估计
summarize(EstMdl);  %展示滞后1阶向量自回归模型模型结果

isStable1=isstable(EstMdl)  %模型稳定检验
% 建立VAR模型，进行检验，得h=1，故VAR模型稳定； 注意这里

YF = forecast(EstMdl,1,ld(n,:));  %一步预测

%% 预测表
fprintf('\n---------------------------------------------------------预测表---------------------------------------------------------\n');
fprintf('%3s%8s%8s%8s%8s\n','数据类型','年份','工业部门','交通运输部门','商业部门');
    fprintf('%3s%14.0f%11.4f%11.4f%14.4f\n','真实值',date(n),ld(n,:));
    fprintf('%3s%14.0f%11.4f%11.4f%14.4f\n\n','预测值',date(n),YF);
